Correlated Gaussian Hyperspherical Method for Few-Body Systems 
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We develop an innovative numerical technique to describe few-body systems. Correlated Gaussian 
basis functions are used to expand the channel functions in the hyperspherical representation. The 
method is proven to be robust and efficient compared to other numerical techniques. The method 
is applied to few-body systems with short range interactions, including several examples for three- 
and four-body systems. Specifically, for the two-component, four-fermion system, we extract the 
coefficients that characterize its behavior at unitarity. 

PACS numbers: 



I. INTRODUCTION 

Ultracold gases in traps or optical lattices have opened 
new possibilities in the study of strongly correlated quan- 
tum systems. From the rich few-body physics of the 
Efimov effect [H 0, H, 0| to the fascinating many-body 
physics of the BCS-BEC crossover [1 i, 0,11, ©, M, El 
E2, Ell E3i experimentalists are now able to realize a 
wide variety of physical systems of great interest for the 
atomic, nuclear, and condensed matter communities. In 
particular, the pureness and controllability of cold atoms 
in optical lattices [HI, Ell E3 make them perfect candi- 
dates for the experimental implementation of condensed 
matter models [see Ref. |l8| and references therein]. In 
all these systems, the rich physics that governs a few 
interacting atoms is crucial for understanding recent ex- 
periments. 

For that reason, extensive efforts have concentrated 
on the development of an accurate description of few- 
body systems. Encouraging advances have been achieved 
in the last decade in the understanding ultracold three- 
body problem d [H, El, H3] • These studies have demon- 
strated the importance of three-body recombination and 
relaxation processes and have determined the effective 
interaction in atom-dimer collisions. Some of these tech- 
niques were subsequently extended to four-body sys- 
tems [HI [H, [H, [24, 25], in a few applications. However, 
the physics of the four-body problem that are far richer 
and more complicated. Also, it is a very challenging 
numerical problem and for that reason it has remained 
largely unsolved except in very limited regimes. Here, 
we present a novel numerical method to handle few-body 
systems that can be used to efficiently describe four-body 
systems, through a combination of different techniques. 

Even though several techniques have been developed in 
recent decades to p rovide solutions for few-body systems 
[H, S3, H, [H, HI , not many of them have been applied 
to numerically solve the Schrodinger equation for systems 
with more than three particles. Among these methods, 
the correlated Gaussian (CG) technique [H [H, [H, [13, 
HE HE H3 m particular has proven to be capable of de- 
scribing a trapped few-body system with short-range in- 
teractions. Because of the simplicity of the matrix ele- 
ment calculation, the CG method provides an accurate 



description of the ground and excited states up to N = 6 
particles [IE [36, 38]. However, the CG method as previ- 
ously implemented can only describe bound states. For 
this reason, previous studies have focused on t rapped sys- 
tems where all the eigenstates are discrete [H|, HI, HI, SO] ■ 
In fact, the CG method requires a nontrivial extension 
in order to describe the continuum and the rich behavior 
of atomic collisions, such as dissociation, rearrangement, 
and recombination processes. 

The hyperspherical representation, in fact, provides 
an appro pria te framework that can treat the contin- 
uum [M Ell, S S H El- In the adiabatic hyper- 
spherical representation, the Hamiltonian is diagonal- 
ized as a function of the hyperradius R, reducing the 
Schrodinger equation to a set of coupled equations in 
a single variable, with a series of different effective po- 
tentials and couplings. The asymptotic behavior of the 
channel potentials describes different dissociation or frag- 
mentation pathways and provides a suitable framework 
for analyzing collision physics. These solutions can be 
readily combined with scattering methods such as the 
R-matrix approach [IE [13, HU to provide an accurate 
description of the collisional dynamics. However, the 
standard hyperspherical methods expand the hyperan- 
gular channel functions in a B-spline or finite element 
basis set [H, HE HH HH, and the calculations become 
very computationally demanding for N > 3 systems. 

It is therefore natural to combine the scalability of the 
CG method with the advantages of the hyperspherical 
representation. In this article, we present an innovative 
way to achieve this combination, in what we term the cor- 
related Gaussian hyperspherical method (CGHS). This 
method uses CG basis functions to expand the channel 
functions in the hyperspherical representation. We show 
that also in this case, the matrix element evaluation is 
greatly simplified thanks to the simple form of the CG 
basis functions. Furthermore, thanks to the explicit cor- 
relation incorporated in these basis functions, only a rel- 
atively small basis set is needed to achieve convergence 
of the lowest channel functions even in the strongly in- 
teracting regime. 

To illustrate the power of the CGHS method, we 
carry out calculations for N — 3, 4-particle systems in 
the strongly interacting regime. First, we analyze sys- 
tems of three-bosons or three fermions at unitarity, and 
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show that the method recovers results that agree with 
semi- analytical predictions. Then, we consider the two- 
component four-fermion system, in the large and positive 
scattering length regime, and reproduce the lowest poten- 
tial curves from Ref. [24(. The CGHS provides a larger 
number of channels which would allow the calculation 
of scattering events not considered in Ref. [24[ . Finally, 
we focus on the universal behavior of four-fermions at 
unitarity. In this regime, the energies of the trapped 
system are trivially determined by the hyperspherical 
potential curves [3a . [53j . Therefore, we can compare 
our calculations with predictions for the trapped sys- 
tem [H, [Hi], [H, [H, H3] • Our results improve and extend 
these previous predictions, and characterize the 20 low- 
est potential curves for even parity and vanishing orbital 
angular momentum. 

This article continues as follows. First, we review both 
the CG and hyperspherical methods in Sec. [Til ln Sec. lIIH 
we introduce the main idea of the CGHS method, leaving 
some details of the implementation for the Appendix [X] 
Sec. IIVI presents our results for three-body systems and 
for the four-fermion system. Finally, Sec. [V] presents our 
conclusions. 



II. THEORETICAL BACKGROUND 

This section discusses the general problem to be solved 
and reviews the correlated Gaussian method. Sub- 
sec. Ill 131 presents the general formalism of the hyper- 
spherical representation and describes how to numeri- 
cally solve the Schrodinger equation in this representa- 
tion using a correlated Gaussian basis set expansion. 

The methods described in this article solve the time- 
independent Schrodinger equation for a Hamiltonian of 
the form 

H = E S v » + V -*W) + E Mnj)- (i) 

where V ex t is an external trapping potential and Vq is 
the interaction potential. The form of the Hamiltonian 
can be varied depending on the particular problem we 
are considering. In the CG method one will usually con- 
sider a spherically-symmetric harmonic trapping poten- 
tial V ext (r) = imiW 2 rf but in hyperspherical calcula- 
tions we usually consider a free system (V ext (r) = 0). 
We can always include the harmonic trapping potential 
in the final step of the hyperspherical calculation, since 
it is a purely hyperradial potential. Depending on the 
particular problem considered, the interacting particles 
will change. For example, all particles interact with each 
other in identical boson systems but only opposite-spin 
fermions interact in two-component Fermi systems (ex- 
cept in a few problems involving p-wave Fano-Feshbach 
resonances). Also, in many cases, the center-of-mass mo- 
tion decouples from the more interesting internal degrees 
of freedom, and it is preferable to use a set of Jacobi 



coordinates rather than the usual single-particle coordi- 
nates. All such options can be treated using the method 
presented below. 



A. Correlated Gaussian Method 

Different types of Gaussian basis functions have long 
been used in many different areas of physics. In partic- 
ular, the usage of Gaussian basis functions is one of the 
key elements of the success of ab initio calculations in 
quantum chemistry. The idea of using an explicitly cor- 
related Gaussian to solve quantum chemistry problems 
was introduced in 1960 by Boys [12] and Singer [3lj . The 
combination of a Gaussian basis and the stochastical vari- 
ational method (SVM) was first introduced by Kukulin 
and Krasnopol'sky [33j in nuclear physics and was exten- 
sively used by Suzuki and Varga [34|, HH, [H, H3] • These 
methods were also used to treat ultracold many- body 
Bose systems by Sorensen, Fedorov and Jensen [58|. A 
detailed discussion of both the SVM and CG methods can 
be found in a thesis of Sorensen [5^1 and, in particular, 
in the book by Suzuki and Varga [27|. In the following, 
we highlight the main ideas of the CG method. 

Consider a set of coordinate vectors that describe the 
system {xi, ...,xjv}- In this method, the eigenstates are 
expanded in a set of basis functions, 

*(x X) --- ,xjv) =^CU$a(xi,--- , X jv). (2) 

A 

Here A specifies a matrix with a particular set of 
parameters that characterize the basis function. It 
is convenient to introduce the following ket notation, 
^(xi, • • • , xat) = (xi, • • • , xjv|^4). Solution of the time- 
independent Schrodinger equation in this basis set re- 
duces the problem to one of diagonalizing the Hamilto- 
nian matrix: 

HC t = E.OC, (3) 

Here, Ei are the energies of the eigenstates, Ci is a vector 
formed with the coefficients Ca and Ti and O are matrices 
whose elements are Hba = (B\H\A) and Oba = (B\A). 
For a 3D system, the evaluation of these matrix elements 
involves 3A^-dimensional integrations which are in gen- 
eral very expensive to compute. Therefore, the effective- 
ness of the basis set expansion method relies mainly on 
the appropriate selection of the basis functions. As we 
will see, the CG basis functions permit fast evaluation of 
overlap and Hamiltonian matrix elements, and they are 
flexible enough to correctly describe physical states. 

To reduce the dimensionality of the problem we can 
take advantage of its symmetry properties. Since the in- 
teractions considered are spherically symmetric, the total 
angular momentum, L, is a good quantum number. For 
simplicity, we will restrict ourselves to L — solutions. 
This restriction allows us to reduce the Hilbert space by 
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introducing restrictions on the basis functions. In partic- 
ular, if the basis functions only depend on the intcrparti- 
cle distances, then Eq. ((2]) can only describe states with 
zero angular momentum and positive parity (L p = + ). 
Furthermore, we can recognize that the center-of-mass 
motion decouples from the system. In such cases, the 
CG basis functions take the form 



<&a(xi, • • • )Xjv) = iPq(Rcm)S { exp 



N 

]>i=i 

(4) 

where S is a symmetrization operator and is the in- 
terparticle distance between particles i and j. Here, 
i/jq is the ground state of the center-of-mass motion. 
For trapped systems, V'o takes the form, ^o(RcAf) = 
e - R cM / 2 ( a ho) . Because of its simple Gaussian form, i/jq 
can be absorbed in the exponential factor. Thus, in a 
more general way, the basis function can be written in 
terms of a matrix A that characterizes them, 



$a(xi,x 2 , ...,xat) = S jexp(-ix T .Ax) 



f 1 " } 

= S j cxp(-- Ai^-Xj) } » ( 5 ) 

where x = {xi, x 2 , xjv}, and A is a symmetric matrix. 
The matrix elements Ay = Aji can be expressed in terms 



of the 



Because of the simplicity of the basis func- 



tions, Eq. (j?]), the matrix elements of the Hamiltonian 
can be calculated analytically. 

The analytical evaluation of the matrix elements is en- 
abled by selecting the set of coordinates that simplifies 
the evaluations. For basis functions of the form of Eq. ([5]), 
the matrix elements are characterized by a matrix M 
in the exponential. Then the matrix element integrand 
greatly simplifies if we rewrite it in terms of the coor- 
dinate vectors that diagonalize that matrix M. This 
change of coordinates permits, in many cases, the an- 
alytical evaluation of the matrix elements. The explicit 
evaluation of several matrix elements can be found in 
Refs. [13, HI]. 

Two properties of the CG method deserve mention at 
this point. First, the CG method does not rely on any 
approximation other than basis set truncation, and the 
solutions can be systematically improved. The accuracy 
of the results are only limited by numerical issues related 
to linear dependence of the basis set. Secondly, the basis 
functions $a are square-integrable only if the matrix A 
is positive definite. This ensures that the wave function 
decays in all degrees of freedom. We can further restrict 
the basis functions by introducing real widths dij such 
that — With this transformation, we ensure 

that A is positive definite. Furthermore, each such width 
is proportional to the mean interparticle distances cov- 
ered by that basis function. Thus, it is relatively easy 



to select the widths after considering the physical length 
scales relevant to the problem. Even though we have re- 
stricted the Hilbert space with this transformation, we 
have numerical evidence that the results converge to the 
exact eigenvalues. 

The linear dependence in the basis set causes prob- 
lems in the numerical diagonalization of the Hamiltonian 
matrix, Eq. To minimize these linear dependence 

problems we restrict the basis function so that the over- 
lap between any two normalized basis functions is below 
some cutoff value. The other method we use to elim- 
inate linear dependence applies a linear transformation 
to produce a smaller orthonormal basis set. 

Finally, we stress the importance of making an appro- 
priate selection of the interaction potential. For the prob- 
lems considered in this article, the interactions are ex- 
pected to be characterized only by the scattering length, 
i.e., to be independent of the shape of the potential. For 
that reason, we can select a model potential that permits 
rapid evaluation of the matrix elements. We have found 
that a model potential with a Gaussian form, 



(0) 



V (r) = -Vbexp ( 



is particularly suitable for this basis set expansion since 
it can be absorbed in the exponential form of the wave 
functions for matrix element evaluation. If the range ro is 
much smaller than the scattering length, then the interac- 
tions are effectively characterized only by the scattering 
length. The scattering length is tuned by changing the 
strength of the interaction potential, Vb, while the range, 
ro, of the interaction potential remains unchanged. This 
is particularly convenient in this method since it implies 
that we only need to evaluate the matrix elements once 
and we can use them to solve the Schrodinger equation 
at any given potential strength (or scattering length). 
Of course, this procedure will give accurate results only 
if the basis set is sufficiently flexible and complete to de- 
scribe the different configurations that appear at different 
scattering lengths. 

In general, this method includes five basic steps: gen- 
eration of the basis set, evaluation of the matrix ele- 
ments, elimination of linear dependence, evaluation of the 
eigenvalue spectrum, followed by a study of stability and 
convergence. The stochastic variational method (SVM) 
Refs. [13, combines the first three of these steps in 
an optimization procedure where the basis functions are 
selected randomly. 



B. Hyperspherical representation 

The main objective of the hyperspherical method is 
to solve the time-independent Schrodinger equation in a 
convenient and efficient way that also provides insight 
into the relevant reaction pathways by which various col- 
lision processes can occur, the first step involves cal- 
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culation of eigenvalues and eigenfunctions of the fixed- 
hyperradius Hamiltonian, which defines the adiabatic 
hyperspherical representation. These eigenvalues and 
eigenfunctions are then used to construct a set of one- 
dimensional coupled equations in the hyperradius R. The 
hyperradius is a collective coordinate related to the to- 
tal moment of inertia of the system [3, In a system 
described by N coordinate vectors ri , . . . , r^v , the hyper- 
radius R is defined by 



[iR 1 



N 
i=l 



(7) 



Here, \x is an arbitrary mass factor called the hyperradial 
reduced mass 6l|and mt are the masses corresponding to 
the particle i. The remaining coordinates are described 
by a set of hyperangles, collectively denoted SI. 

The total number of spatial dimensions of this N- 
particle system is d = 3 N. The total wave function 
tp is rescaled by R, $ = R^ d ~ l 'i 2 ij), so that the hy- 
perradial equation resembles a coupled one-dimensional 
Schrodinger equation. In the adiabatic representation, 
the wave function ^ E (R, Si) is expanded in terms of a 
complete orthonormal set of angular wave functions 
and radial wave functions F vE , such that 



v E (R, si) = X Fue{R)$»{R; si). 



(8) 



The adiabatic eigenfunctions, or channel functions <!>„, 
depend parametrically on R and are eigenfunctions of a 
3 iV — 1 partial differential equation (which reduces to 
3iV — 4 dimensions if the center-of-mass motion is re- 
moved explicitly): 



2A2 



K 2 K 
\2[iR 2 



{d-l){d-3)h 2 
8fiR 2 



V(R,Sl)j $ v (R;(l) 

= U v {R)$ u {R;Sl). (9) 



Here, A is the grand angular momentum operator, which 
is related to the kinetic term by 



E 



2m, 



h 2 1 d d _ x d A 2 h 2 
2^R d ^m R dR + 2^R 2 ' 



(10) 



The U U {R) obtained in Eq. §§§ are effective hyperra- 
dial potential curves that appear in a set of coupled one- 
dimensional differential equations: 



h 2 d 2 



2/x dR 2 



U„{R) 



F uE (R) 



2P^(R)— + Q^(R) 
dR 



F u , E {R)=EF uE {R). 

(11) 

These differential equations [Eq. (fTTj) ] are coupled 



through the P vv i (R) and Q vv i (R) couplings defined as 



P UV ,(R) = (^(R;Sl)\-^\^ v ,{R;Sl)) 
Q VV .{R) = ($ u (R;Sl)\-^\$ v ,(R;Sl)) 



(12) 
(13) 



Since the basis set expansion of the wave function, 
Eq. ([8|), is complete in the 3 iV-dimensional space, 
Eqs. ^ and (fTTj) reproduce exactly the original d- 
dimensional Schrodinger equation. As in most numeri- 
cal methods, the solutions are approximated by truncat- 
ing the Hilbert space. In this case, the Hilbert space is 
truncated by considering a finite number of channels in 
Eq. (|11[) . This approximation is easily tested by analyz- 
ing convergence with respect to the number of channels 
included in the calculation. 

The utility of the hyperspherical representation relies 
on the assumption that the wavefunction variation with 
the hyperradius R is smooth. In such cases, only a few 
channels are relevant, and the couplings are small and 
vary smoothly with R. Furthermore, a fairly good ap- 
proximation to the solutions can be achieved by truncat- 
ing the expansion in Eq. (jSj) to a single term: 



* E {R,Si) = F vE {R)$ v {R-Si). 



(14) 



This adiabatic hyperspherical approximation leads to an 
effective one-dimensional Schrodinger equation, 



h 2 d 2 

2[i dR 2 



W V {R) 



F vE (R) — EF vE (R) , 



where the effective potential is 



W V {R) = U V {R) - —Q VV {R). 

2/i 



(15) 



(16) 



Here, the first term is the hyperradial potential curve, 
and the second term is "adiabatic correction" , i.e. the 
repulsive kinetic contribution of the hyperradial depen- 
dence of the channel function. If the potential curves are 
well-separated and have no strong avoided crossings in 
the relevant range of energy and radius, then the adia- 
batic approximation can be quite accurate for the lower 
states in any given potential curve. This approximation 
comes from a truncation of the Hilbert space and, for 
that reason, obeys the variational principle. Any dis- 
crete energy eigenvalue obtained with this method is an 
upper bound of the exact energy level, in the sense of the 
Hylleraas-Undheim theorem. An approximate descrip- 
tion of the spectrum can be achieved by combining the 
energies obtained from the adiabatic approximation ap- 
plied to each channel separately. For example, bound 
states of excited potential curves which are above the 
lowest fragmentation threshold would represent quasi- 
bound states. This approach is equivalent to neglecting 
all off-diagonal couplings in Eq. (fill) an d produces an 
approximate spectrum which is not variational. Another 
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useful approximation is obtained by neglecting the sec- 
ond term in Eq. (fTB]) , i.e., replacing W (R) by U (R) in 
Eq. (fT5| . This is usually called the hyperspherical Born- 
Oppenhcimer approximation. As in the standard Born- 
Oppenheimer approximation for a diatomic molecule, the 
approximate energy obtained in this manner re pre sents 
a lower bound to the exact ground state energy [621 ]. 

Next, we show how Eq. © is solved, and how the P vv i 
and Q vv i are evaluated. 



where 



(24) 



Thus, we can obtain all the couplings from the evalua- 
tion of P and Q. In the basis set expansion, P and Q 
can be calculated using matrix multiplication. With the 
expansion in Eq. (|19p. 



(25) 



\%(R)) = '£\B i )c ip + \B i )c il 



C. Expansion of the channel function in a basis set 



In the hyperspherical method (see Sec. IIIBj ), channel 
functions are eigenfunctions of the adiabatic Hamiltonian 

H A (R; si), 



H A (R; tt)$v(R; SI) = U V (R)$ V {R; SI). 



(17) 



The eigenvalues of this equation are the hyperspherical 
potential curves U V (R), which serve as readily visualiz- 
able reaction pathways. The adiabatic Hamiltonian has 
the form: 



H A (R;Sl) = 



h 2 A 2 (d-i){d-3)h 2 

2^iR 2 + 8fiR 2 



V(R,Sl). (18) 



Here, d = 3Nj where Nj is the number of Jacobi coor- 
dinate vectors. 

A standard way to solve Eq. (fT7j) is to expand the 
channel functions in a basis, 



\$»(R;Sl)) =J2\B t (R;Sl))^(R). 



(19) 



Here /z labels the channel function. The \Bi(R;Sl)) are 
the basis functions. With this expansion, Eq. (TTTj) re- 
duces to the eigenvalue equation 



H A (RK = U^(R)0(R)c, 



(20) 



The [i-th column vector = {c^j-ji = 1, ...£), where 
D is the dimension of the basis set. Ha and O are the 
Hamiltonian and overlap matrices whose matrix elements 
are given by 



n A (R)ij = (BilHAiR-^Bj) 
0{R) l3 = (Bi\Bj) 



(21) 
(22) 



Once the hyperradial potential curves are calculated, 
we still need to evaluate the P and Q non-adiabatic 
couplings between the channel functions (Eq. (Tj"2")) in 
Sec. Ill Bp . To evaluate the Q coupling, we use the iden- 
tity 



Qv^{R) = —Qv^{R) 



dP^(R) 
OR ' 



(23) 



Here and in the following, we have omitted the radial and 
angular dependence of functions, and we have introduced 
the notation F for the derivative of F with respect to R. 
The P coupling takes the form 

Pun = E <Zi ( B i\ B J <V +C ^ ( B il^*> c ^ = clO^+clVc^. 

(26) 

where V{R) is defined later in Eq. (|2"9")) . The same pro- 
cedure can be done for the Q matrix elements with 



+ c T V) (B^) 6i M + c T VJ (Bj\Bi) c v (27) 

and can also be written in terms of matrix multiplica- 
tions: 

- ?lO{R% + ?lV{R)c^ + clV T {R% + clQ{R)^. 

(28) 

In Eqs. (|26l |2"8|) we have used the overlap matrix O and 
defined the matrices V and Q whose matrix elements are 

V{R)a = (Bi(R)\Bj) and Q(R) l3 = (B^) . (29) 

The derivatives of the (R) coefficients that form the 
are calculated numerically using the three point rule. 



III. CORRELATED GAUSSIAN 
HYPERSPHERICAL METHOD 

As we have seen in the previous section, the implemen- 
tation of hyperspherical calculations requires the eval- 
uation of the Hamiltonian matrix elements at fixed R 
(Eqs. |2"T1 and |2"2"|) . This is one of most time consuming 
part of the calculation which for an N — 4 system re- 
quires a 5 dimensional integration. Thus, we need to 
find an efficient way to evaluate Hamiltonian matrix el- 
ements at fixed R. As a prelude, we first review how 
multidimensional matrix elements evaluations reduce to 
analytical forms in the standard CG method. This will 
be the key to evaluating matrix elements in the hyper- 
spherical variant of this method. 

In the CG method, we select, for each matrix element 
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evaluation, a set of coordinate vectors that simplifies the 
integration, i.e., the set of coordinate vectors that diago- 
nalize the basis matrix M which characterizes the matrix 
element. The flexibility to choose the best set of coordi- 
nate vectors for each matrix element evaluation is crucial 
for the economy of the CG method. 



This selection of the optimal set of coordinate vectors 
is formally applied by an orthogonal transformation from 
an initial set of vectors x = {xi,...,xn} to a final set 
of vectors y = {yi, yjv}: Tx = y, where T is the 
orthogonal transformation matrix. The hyperspherical 
method is particularly suitable for such orthogonal trans- 
formations because the hyperradius R is an invariant un- 
der them. Consider the hyperradius defined in terms 
of a set of mass-scaled Jacobi vectors 0, H3, [5^, l63| . 
x = {xi, ...,x N }, 



fxR 2 = /x^x; 2 , 



(30) 



If we applied an orthogonal transformation to a new set 
of vectors y, then 



fiR 2 = fi ^2 x i 2 = w TT Ty = a* ^2 yi 2 



(31) 



where we have used the fact that T T T = 7, and I is the 
identity. Therefore, in the hyperspherical framework we 
can also select the most convenient set of coordinate vec- 
tors for each matrix element evaluation. This will be the 
key to reducing the dimensionality of the matrix element 
integration. This transformation amounts to selecting, 
for each matrix element evaluation, the set of hyperan- 
gles (fi) that simplifies the matrix-element evaluation. 



As an example of how the dimensionality of matrix- 
element integration is thereby reduced, consider an L = 
three-dimensional A^-particle system with the center of 
mass removed. It can be shown that this technique re- 
duces a (3N — 7) numerical integration [H| to a sum 
over the symmetrization permutation of (N — 3) numer- 
ical integrations. This result implies that for N = 3 the 
matrix element evaluation can be done analytically (see 
Appendix IA 1[) and that for N = 4, it requires a sum of 
one-dimensional numerical integrations [651 ]. 



Once the basic idea of the appropriate change of vari- 
ables for each matrix element calculation is understood, 
the actual calculation of the matrix elements using cor- 
related Gaussian basis function is straightforward. Ap- 
pendix IA II shows, as an example, how the matrix ele- 
ments can be calculated analytically for a three particle 
system (the calculation of the matrix elements for N = 4 
are not presented here but can be found in Ref. [65[). Fi- 
nally, Appendix IA 21 discuss in general how this method 
is implemented. 



IV. RESULTS 

In this section, we present CGHS results for iV = 3, 4. 
First, we analyze two different N = 3 systems and com- 
pare them with analytical predictions. Then, we present 
four-fermion potential curves and compare them with re- 
cent predictions [24J. Finally, we characterize the four- 
fermion L = potential curves at unitarity and extract 
the s„ coefficient that characterize the universal regime. 

To test the CGHS method, we calculate the hyper- 
spherical potential curves at unitarity for three interact- 
ing bosons. For zero-range interactions, the potential 
curves at unitarity are inversely proportional to the hy- 
perradius. For example, the lowest potential curve for 
three identical bosons is given by 



U (R) = -'- 



1/4 



2fiR 2 



(32) 



The coefficient sq sa 1.0062 can be obtained analytically 
in the theory of Efimov states [HI3.[66|. A simple and fast 
numerical CGHS calculation with only 30 basis functions 
extended up to R — lOOro shows, at large R, the expected 
1/R? behavior. Extrapolation of our potential curves to 
R — > oo gives so « 1.0059. 

Similarly, we analyze the system of two indistinguish- 
able fermions resonantly interacting with a third particle 
of equal mass. For such system, the zero-range model 
predicts a lowest potential of the form, 



Uo(R) 



Pi 1/4 
2fiR 2 



(33) 



The value of po w 2.166222 can also be predicted analyti- 
cally. Using a slightly larger basis set of 90 basis function 
we extend the CGHS calculations up to R — 4000ro. Ex- 
trapolating our potential curves to R — > oo we obtain 
po « 2.166218. 

These two examples show that the CGHS method is 
flexible enough to describe a strongly interacting system 
with relatively small basis sets and analytical matrix el- 
ement evaluations. The main limitation of these calcula- 
tions come from linear dependence issues. At the A^ = 3 
level, this method cannot probably compete with more 
sophisticated calculations which permit calculations up 
to R = 10 6 ro [3, 52]. However, it has been a challenge to 
extend hyperspherical methods beyond A^ = 3. One suc- 
cessful method uses Monte Carlo techniques to describe 
the lowest channel function and extends it application to 
large (N < 10) systems [67|. However, this method can 
only calculate the lowest potential curve and leads to an 
approximate solution. In contrast, the CGHS method 
can be readily extended to A^ = 4 particles (and possi- 
bly beyond) and allows to obtain a full solution which 
represents the current state of the art of hyperspherical 
methods. 

The development of four-body hyperspherical methods 
allows, for one thing, an analysis of the full energy de- 
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TABLE I: Non interacting coefficients p v of the four-fermion 
potential curves and their degeneracies A„. 



400 600 800 

R/r 



1000 1200 



FIG. 1: (Color online) Adiabatic hyperspherical potential 
curves U V (R) (solid lines) for two spin-up and two spin-down 
fermions with an atom-atom scattering length a s — 100ro. 
The dashed line at E — 2Et (blue online) is the dimer-dimer 
threshold, the dashed line at E = Et (red online) is the dimer- 
two-atom threshold, and the dashed line at E = (green 
online) is the four-atom threshold. Dashed curves are predic- 
tions from Ref. 0. 



pendence of the dimer-dimer scattering length. Figure [T] 
presents the four-fermion potential curves obtained with 
the correlated Gaussian hyperspherical mcthod(CGHS). 
There are three relevant energy thresholds mark with 
dashed lines in Fig. [T] dimer-dimer threshold at 2E],, 
dimer-two-atom threshold at E\, and four-atom thresh- 
old at energy. The lowest curve represents the dimer- 
dimer channel and potential curves going asymptotically 
to Eb and represent dimer-two-atom and four-atom 
channels, respectively. Standard multichannel scatter- 
ing techniques, like the R-matrix method, can be applied 
to solve the hyperspherical coupled differential equations. 
This analysis was performed in a recent study by D'Incao 
et. al. [2J| , which obtained the energy dependence of the 
dimer-dimer scattering length for equal mass systems. 
Black dashed curves in Figure Q] represent the potential 
curves of Ref. (24[. As we can see, the CGHS method pre- 
sented here predicts very similar potential curves. The 
dimer-dimer potential curves obtained with the differ- 
ent methods are almost indistinguishable. For dimer- 
two-atom potential curves, the CGHS predicts lower po- 
tential curves suggesting that the CGHS calculation is 
slightly better. At large R, the asymptotic behavior of 
both methods agree. This is very encouraging since in 
the method of D'Incao et. al., the asymptotic behavior of 
the channel functions is correct by construction, whereas 
in the CGHS it constitutes an important, nontrivial 
test. Preliminary calculations with the CGHS potential 
curves predict a similar energy dependence of the dimer- 
dimer scattering length. Therefore, the CGHS opens the 
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possibility for accurately analyzing four-body scattering 
events, as has been carried out for four-interacting bosons 
in Refs. [H[6|. 

The calculations of the potential curves at unitarity al- 
lows us to extract the four-fermion universal coefficients. 
As in the N = 3 system, the potential curves can be 
written as 153, [701 , 



UAR) = 



1/4 



2fiR 2 



(34) 



This functional form of the potential curves was verified 
indirectly in Ref. [38[ by analyzing the spectrum of the 
four-fermion system under spherical harmonic confine- 
ment. It can be shown that all the couplings vanish when 
the potential curves are proportional to 1/R 2 . There- 
fore the system is described by a set of uncoupled one- 
dimensional Schrodinger equations that can be solved an- 
alytically once the trapping potential is included. These 
procedure leads to simple expressions for the trapped en- 
ergies [53| 



E vn = {Pu + 2n + 5/2)Hu, 



(35) 



where u> is the trapping frequency and we have included 
the zero point energy of the center-of-mass motion. In 
Ref. the 2hw spacing was verified and the low- 

est p v coefficients were identified. Equations (f3~4"|) and 
(|35[) are also valid in the non- interacting limit. For 
L = and positive parity solutions, the p^ 1 values and 
their degeneracies \ v have relatively simple closed forms: 
pff 1 = 11/2 + 2v and A„ = i/ 4 /96 + 7^ 3 /48 + 17^ 2 /24 + 
133^/96 + 57/64 + (-l) v v/32 + 7/64(-l) ,y . Their lowest 
values can be found in Table |TJ 

The development of the CGHS method allows us 
to carry out a hyperspherical calculation for the four- 
fermion problem and to directly verify the form of the 
hyperspherical potentials. Also, it allows us to analyze 
deviations from the zero-range solutions due to finite- 
range effects. 

The 20 lowest four-body potential curves 
2iiR 2 U v {R)/h 2 for the equal-mass system are pre- 
sented in Figure [U We can identify three regimes 
in these potential curves. The region R < r$ is 
controlled by the kinetic energy. The kinetic energy 
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10' 

R/r 



FIG. 2: (Color online) Hyperspherical potential curves at 
unitarity (a — > oo) for the 4-fermion system multiplied by 
2fj,R 2 /h 2 . The solid lines represent the predictions from ana- 
lyzing the spectrum obtained with the CG method. The sym- 
bols correspond to direct evaluation of the potential curves 
with the CGHS method. 



effects are more important than the interaction energy 
and the potential curves are well approximated by 
the non-interacting potential curves. In other words, 
2/iR 2 U v (R)/h 2 fn (p^ 1 ) 2 - 1/4 and the eigenchannels are 
well approximated by the hyperspherical harmonics (see 
Sec. Ill B| ). For that reason, there is a large degeneracy in 
the R < tq region which corresponds to the degeneracy 
of the A 2 operator. Furthermore, the potential curves 
are, to a good approximation, proportional to 1/R 2 . 
The second region is ro < R < 20ro. In this region 
both the kinetic and the interaction terms are important 
and finite range effects are important. In the third 
region, R > 20ro, the potential curves recover their 
universal behavior. The potential curves are, again, 
approximately proportional to 1/R 2 . As R/ra increases, 
finite-range effects tend to zero and we obtain the 
zero-range potential curves at unitarity. Therefore, 
in this region, the eigenvalues of 2nR 2 U u (R) /ti 2 are 
approximately (p 2 — 1/4). Thus, we can compare these 
results with the ones deduced from trapped calculations 
for ro/aho = 0.01 presented in Ref. [38|. The solid lines 
correspond to (p 2 - 1/4), (p\ - 1/4) and (p 2 - 1/4) 
respectively 71[. There is good agreement between 
the predictions from the trapped system obtained with 
CG and the direct computation of the potential curves 
through CGHS. 



To quantify this last statement, we analyze the value 
of p . Several groups [H, 0, [H, HE H3 have tried to 
benchmark the four-body value of Eqq, which is simply 
related to po. The calculations from Ref. [5(| use zero- 
range interactions explicitly and they report a value of 
E 00 w (5.045 ± 0.003)huj. To extract the p value in 



the zero-range limit, we carry out two different calcula- 
tions. First, we study the Eqq energy obtained with the 
standard CG method as a function of the range of the 
two-body interaction and then we extrapolate to zero- 
range limit. This method was previously applied for 
the three-body system and the numerical results agreed 
with the analytical predictions up to 7 digits [40]. The 
same procedure applied to the four-body system, leads 
to po w 2.5096. The second calculation analyzes the 
long-range behavior of the potential curves. To eliminate 
finite-range effects, we extrapolate the potential curve 
Uo(R) to R/ro — ► oo. In this limit, Uq(R) is character- 
ized by a value po « 2.5092. These two different methods 
provide a value of po which agrees in four digits. These 
values are slightly lower than po ss 2.545 ± 0.003 pre- 
dicted in Ref. (5(| • This suggests that the uncertainty in 
Ref. [56[ was apparently underestimated. 

The calculations of the lowest 20 universal coefficients 
p v is reported in Table [Til It is interesting to note that 
some of the p v coefficients are very similar to the non- 
interacting coefficients. For example, the p v coefficient 
for u=12, 13, 14 coincide with noninteracting p^ 1 coef- 
ficient. Two of these potential curves are also described 
by Pv — 9.5 in the small R region and deviate from these 
value in the region R ~ ro . These channels have nodes in 
every spin-up-spin-down interparticle distance, therefore 
at large distances they recover the noninteracting behav- 
ior. The third potential curve smoothly decrease from 
pfj 1 = 11.5 at small R to p v = 9.5 at large R. 

Finally, note that the CGHS method has been suc- 
cessfully applied to the four-boson system [68|. In that 
study, the four-boson spectrum is calculated from the 
CGHS potential curves. Also, that study considers scat- 
tering events such as four-body recombination, which was 
calculated and predicted to be important for the under- 
standing of a recent experiment on Efimov physics in an 
ultracold Bose gas 
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V. CONCLUSIONS 

We have presented an innovative numerical method 
suitable for the analysis of four-body processes. We 
have shown several examples for three- and four-particle 
systems recovering known results. Furthermore, we 
have obtained the lowest 20 p v coefficients for the two- 
component four-fermion system that characterize both 
free and trapped systems. These coefficients also charac- 
terize the spectrum of four trapped fermions at unitar- 
ity. Our results considerably extend previous calculations 
and provide more accurate energies. 

The CGHS method has been used to analyze the four- 
boson system in Refs. [H, [H, I72 I . p redicting new phe- 
nomena observed experimentally [73j ■ It has also built a 
theoretical foundation for the analysis of four-body colli- 
sional processes in other s yste ms such as two-component 
Bos e sy stems [7J, [75|, [zl, |77j, Bose-Fermi mixtures l78l. 
l79l . Ho], and three-component Fermi gases [8l|, HH, |83| |. 
Even though this method was initially implemented to 
treat ultracold systems using model potentials, it can be 
in principle extended to other four-body problems. 

The authors would like to thank S. T. Rittenhouse, N. 
P. Mehta and J. P. D'Incao for useful discussions and 
for providing their four-fermion numerical data (dashed 
curves included in Fig. [I]). This work was supported in 
part by NSF. 



For equal mass systems, we can write Eq. (|All) in terms 
of the following Jacobi coordinates: 



xi 



V2 



(n - r 2 ), 



x 2 



r.3 



ri + r 2 



3 V 2 
The basis functions [Eq. (|A1[) ] can be written as 



(A2) 
(A3) 



*A(ri2,ri3,r 23 ) = (xi,x 2 |A) = exp(- 



x T i.x, 



exp(- 



Xi.Xidn + 2xi.x 2 ai2 + x 2 .x 2 a 22 



) (A4) 



where x = {xi,x 2 } and A is a 2 by 2 symmetric matrix 
whose elements are an = 2/rff 2 + 1/2(1/^13 + V^ls) > 
a X2 = oai = \/3/2(l/dl 3 -l/df 3 ), and a 22 = 3/2(l/d? 3 + 
l/<i 23 ). In Eq. (|A4I) . we can clearly see that the state 
(xi, x 2 |A) depends only on the distances x±, and x 2 plus 
the angle #i 2 between them, cos#i 2 = Xi.x 2 /a;ia; 2 . 

We want to obtain the matrix elements corresponding 
to these basis functions at fixed hyperradius R. We define 
the hyperradius to be R 2 — x. 2 + x\ . The integrand of 
the overlap matrix element between \A) and \B), noted 
as B.A, is 



B.A 



exp( 



x T .(A + B).x 



(A5) 



APPENDIX A: APPLICATION OF 
CORRELATED GAUSSIANS TO THE 
HYPERSPHERICAL FRAMEWORK 



This appendix illustrates how correlated Gaussian ba- 
sis functions can be used in the general hyperspherical 
framework presented in Subsec. Ill CI First, we consider 
the three-particle case and calculate the matrix elements 
fEqs. |2"T1 12"21 and |2"9"]) . Then, we discuss how to generate 
and optimize the basis set. 



We change to the Jacobi basis set that diagonalizes A+B 
and we call /3i and /3 2 the eigenvalues and y = {yi,y 2 } 
the orthonormal eigenvectors. In this new coordinate 
basis, Eq. (|A5|) has a simple form, 



B.A = exp(- 



)■ 



(A6) 



We integrate over the angles of the vectors yi and y 2 and 
we fix the hyperradius, so y\ = Rcos8 and y 2 = RsmO. 
In this set of coordinates, the matrix element at fixed R 



is 



Unsymmetrized matrix element evaluation for 
three particles 



In this subsection, we present as an example the eval- 
uation of the matrix elements (Eqs. [21] [22j [29]) of three 
particle system. Consider a system in which the center- 
of-mass motion decouples. Then, the L p = + solutions 
of the body-fixed system can be expanded in terms of 
the interparticle distances. For the three-body system 
the correlated basis functions take the form 



^A(ri2,r 13 ,r 23 ) = exp 



' i 2 
2d 2 2 



1 13 

2d 13 



' 23 
2d 2 3 



(Al) 



B\A) = {At:) 2 



tt/2 



-R (01 cos 2 8+f} 2 sin 2 9)/2 



COS 



II 



This integration has a closed- form result, 

3 exp(-ft±^ 

= 2tt 

R £ 



(B\A) 



MO 



(A7) 



(A8) 



Here we have introduced the definition £ = R 2 (f3\ — /3 2 )/4. 

To simplify the interaction matrix element evaluation, 
we can adopt a Gaussian model potential as was utilized 
in the CG method. In this case, the interaction term 
can be evaluated in the same way we have calculated 
the overlap term since the interaction is also a Gaussian. 
Each pairwise interaction can be easily written as Vij = 
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V exp( 



2r?, 



y exp(— x T .Af( i:, '.x/(2r 2 )). Therefore, to the basis set that diagonalizes A + B. We obtain 



to calculate the interaction matrix element, we need to 
evaluate 



(B\V ij \A)=V / drtexp( 



.(A + B + MW/r%).x. 



(A9) 

This integration can be done following the same steps 
of the overlap matrix element. Equation (|A8|) can be 
used directly if we multiply it by Vo, and fa and fa are 
replaced by the eigenvalues of A + B + jr\. Note 
that for each pairwise interaction (and for each pair of 
basis functions in the matrix element), the matrix 
changes and requires a new evaluation of the eigenvalues. 

The third term we need to evaluate is the hyperangular 
kinetic term at fixed R. This kinetic term is proportional 
to the grand angular momentum operator A defined for 
the N = 3 case as 



A 2 h 2 
2^R 2 



E 



2^ 



h 2 1 d * d . 



The expression can be formally written as 



To, =Tt — Tr 



where 



and 



To = 



A 2 h 2 
2{iR 2 ' 



T T = - 



v-£ 2 V 2 



2^ ' 



2^R 5 dR dR' 



(A12) 



(A13) 



In typical calculations, Tq is evaluated by directly ap- 
plying the corresponding derivatives in the hyperangles 
fi. However, in this case, it is convenient to evaluate 7~t 
and Tr separately, and make use of (|A11[) . 

The integrand of the total kinetic term Tt takes the 
form 



B\T T \A = exp(- 



x 1 .B.x 



) -E 



2[L 



exp(- 



x .Ax 



2 

(A14) 

First, we diagonalize A and use the eigenvectors and 
eigenvalues of A, obtaining, 

B\T T \A = (-Tr[A] + x T .A 2 .x) exp( ^ T(i + B) ~ 



)• 



(A15) 

Here Tr is the trace function. We can use TV [A] = (ct\ + 
02), where a\ and 012 are the eigenvalues of A. Now 
we diagonalize A + B. We call T the matrix with the 
orthonormal eigenstates in columns and fa and fa & r e the 
eigenvalues of A + B. We make a change of coordinates 



B\T T \A = -g (-3K + a a ) + y.G.y) exp(-M±M 

^ (A16) 
where G = T T .A 2 .T, and yi and y2 are the vectors in 
the new eigen basis. The integration over the angles of 
these vectors is trivial. After this integration, we fix the 
hyperradius and integrate over the hyperangle 9 defined 
by Ui = RcosO and 1/2 — Rsin9, 



(B\T T \A) 



(4ir) 2 h 2 r' 2 r 



2/y 



r r 

Jo L 



gnR 2 cos 2 9 + g 2 2R 2 sin 2 9 



exp(- 



faR 2 cos 2 9 + faR 2 sin 2 9 



) cos 2 9 sin 2 0d0. (A17) 



This integration can be done analytically and the results 
expressed in terms of the Bessel functions I\ and Iq: 



(B\T T \A) 



„ Oi+fei-R 2 3d2 



( AU ) |-8( 5 ii ~ .92 2 )/o K] + |{8(ffii - 522) 



+(/3i - /3 2 ) (-6(ai + a a ) + On + S 22 )i? 2 ) }/i [£]} • 

(A18) 

Now we will evaluate 7r, the hyperradial kinetic term. 
It is written as 



fi 2fi\R 5 / 2 dR 2 AR 2 
Therefore, the integrand takes the form 



(A19) 



B\T R \A = - 
( 1 d 2 



h 2 , x T .B.x, 
^exp( — ) 



\R 5 / 2 dR 2 



AR 2 ) 



R 5 ^ 2 I <'xp 



x .Ax 



(A20) 



We use the property x T .Ax = R 2 Fa(Q) to evaluate 
the derivatives with respect to R. This allows a simple 
calculation of the derivatives in Eq. (|A20p , yielding 

B \T R \A = ^ [6x T .Ax- (x T .Ax) 2 ] e~^- ( A + B )-*/ 2 . 

(A21) 

Next we diagonalize A + B, and set D = T T .A.T, giving 

B\T R \A = —2 L6y.Z?.y - (y.D.y) J exp I 

(A22) 

The terms y.D.y and (y.D.y) 2 depend on the polar an- 
gles of the vectors. The integration over the polar angles 
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(ill = {01, ^1} and il 2 = {<fe, #2}) of these terms is 



General considerations 



J [Zy.D.y-{y.D.yf] drifts = (4^) 2 {&d llV l+M 22 y 2 2 

- [d 2 n yf + (2d u d 22 + 4d 2 2 /3)2/ 2 t/ 2 + d^vi] }• (A23) 

Now we carry out the integration over the hyperangle 9, 
using yi = RcosO and y 2 = i?sin0, which gives 

+ 6d 22 i? 2 sin 2 6-d^R 4 cos 4 

- (2d u d 22 + 4d 2 2 /3)i? 4 cos 2 6> sin 2 - dj 2 R 4 sin 4 0} 
ftj? 2 cos 2 + /3 2 i? 2 sin 2 



exp 



cos 2 0sin^ 0d0. 

(A24) 



This integration has the analytical form 



(B\T R \A) 



4 7T R 



R fl 64£ 2 I 

+ (iii - i 22 )(6(-A + fo + in ~ ^22)+ 



8df 



4?(d 11 + d 22 ))J/ [e] + ^ 



64d 2 5 



48(dn - d 22 )(-Pi +f3 2 + in - d 22 ) 
+ 8f(-3/9i + 3/3 2 + 4d u - 4i 22 )(in + d 22 ) 

+16^ 2 (d 2 1 +d 2 2 )]/ 1 K]}. (A25) 



Combining Eqs. (|A181 IA25|) . we obtain 7h. The 
expression for To can be simplified using the relation 
G = D 2 to write G matrix elements of Eq. (|A18|) in terms 
of the ones of D. This same procedure can be applied to 
extract the P and Q matrix elements: 



(B 



OA 
'dR' 



and 



'dR^dR,' 



(A26) 



A useful test to verify the functional form of the ma- 
trix elements is to integrate them with respect to R, with 
the corresponding volume clement, and compare that re- 
sult with the standard CG matrix elements. Another 
important test is to verify that Tq is symmetric under 
the exchange of the basis functions A and B. This is not 
a trivial test since neither Tt nor 7r are symmetric. 

A major advantage of these matrix element evaluations 
is that they can be easily extended to four particles. In 
general, these matrix elements evaluations would require 
a 5-D numerical integration but for these basis functions, 
with the above analytical development, they only require 
a TD numerical integration. 



Many of the procedures of the standard CG method 
can be easily extended to the CGHS. The selection, sym- 
metrization, and optimization of a basis follow the same 
ideas of the standard CG method. However, the evalua- 
tion of the unsymmetrizcd matrix elements at fixed R is 
clearly different. Furthermore, the hyperangular Hamil- 
tonian [Eq. [T7] need to solved at different hyperradius 
R. 

There are several properties that make this method 
particularly efficient. For the model potential used, the 
scattering length is tuned by varying the potential depths 
of the two-body interaction. Therefore, as in the CG 
case, the matrix elements need only be calculated once; 
then they can be used for a wide range of scattering 
lengths. Of course, the basis set should be complete 
enough to describe the relevant potential curves at all 
the desired scattering length values. 

The selection of the basis function generally depends 
on R. To avoid numerical problems, the mean hyperra- 
dius of each basis function (R) B should be comparable 
to the hyperradius R in which the matrix elements are 
evaluated. We can ensure that (R) B ~ R by selecting 
some (or all) the weights d%j to be of the order of R. 

We consider two different optimization procedures. 
The first possible optimization procedure is the following: 
First, we select a few basis functions and optimize them 
to describe the lowest hyperspherical harmonics. The 
Gaussian widths of these basis functions are rescaled by 
R at each hyperradius so that they represent the hyper- 
spherical harmonics equally well. These basis functions 
are used at all R, while the remaining arc optimized at 
each R. Starting from small R (of the order of the range 
of the potential), we optimize a set of basis functions. As 
R is increased, the basis set is increased and reoptimized. 
At every R step, only a fraction of the basis set is opti- 
mized, and those basis functions are selected randomly. 
After a several i?-steps, the basis set is increased. 

Instead of optimizing the basis set at each R, one 
can alternatively try to create a complete basis set at 
large R m ax- In this case, the basis functions should be 
complete enough to describe the lowest channel func- 
tions with interparticle distances varying from interac- 
tion range ro up to the hyperradius R m ax- Such a basis 
set can be rescaled to any R < R max and should effi- 
ciently describe the channel functions at that R. The 
rescaling procedure is simply d%j/R = d^ ax / R max . This 
procedure avoids the optimization at each R. Further- 
more, the kinetic, overlap, and couplings matrix elements 
at R are straightforwardly related with the ones at R ma x ■ 
So, the interaction potential is the only matrix element 
that needs to be recalculated at each R. This property 
can be understood using dimensional analysis. The ki- 
netic, overlap, and coupling matrix elements only depend 
on R, so a rescaling of the widths is simply related to a 
rescaling of the matrix elements. In contrast, the interac- 
tion potential introduces a new length scale, so the ma- 
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trix elements depend on both R and r , and the rescaling 
does not work. 

These two methods, the "complete basis set" or the 
"small optimized basis set" method, can be appropriate 
in different circumstances. If a large number of channels 
are needed, probably the complete basis method is the 
best choice. But, if only a couple of particular channel 
potential curves and couplings are needed, then the small 
optimized basis set method might be more efficient. 

The most convenient strategy we have found for op- 
timizing the basis function in the four-boson and four- 
fcrmion problems is the following: First we select an 
hyperradius R m that is R m w 300 tq where the basis 
function will be initially optimized. The basis set is in- 
creased and optimized until the relevant potential curves 
are converged and, in that sense, the basis is complete. 



This basis is then rescaled, as proposed in the second op- 
timization method, to all R < R m . For R > R mi it is 
too expensive to have a "complete" basis set. For that 
reason, we use the "small optimized basis set" method 
which allows a reliable description of the lowest poten- 
tial curves. 

Note that for standard correlated Gaussian calcula- 
tions, the matrices A and B need to be positive definite. 
This condition restricts the Hilbert space to exponen- 
tially decaying functions. In the hypersphcrical treat- 
ment, this is not necessary since the matrix elements 
can always be calculated at fixed R, as the integrals 
converge even for exponentially growing functions. This 
gives more flexibility in choosing the optimal basis func- 
tions. 



[1] E. Braaten and H. W. Hammer, Phys. Rep. 428, 259 
(2006). 

[2] V. Efimov, Yad. Fiz. 12, 1080 (1970 [Sov. J. Nucl. Phys. 

12, 589 (1971)]). 
[3] B. D. Esry, C. H. Greene, and J. P. Burke Jr, Phys. Rev. 

Lett. 83, 1751 (1999). 
[4] T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. 

Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, 

H. C. Naegerl, et al., Nature 440, 315 (2006). 
[5] D. M. Eagles, Phys. Rev. 186, 456 (1969). 

[6] A. J. Leggett, in "Modern Trends in the Theory of 
Condensed Matter", ed. by A. Pekalski J. Przystawa, 
Springer, Berlin (1980). 

[7] P. Nozieres and S. Schmitt-Rink, J. Low Temp. Phys. 59, 
195 (1985). 

[8] M. Greiner, C. A. Regal, and D. S. Jin, Nature 426, 537 

(2003) . 

[9] S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. 
Riedl, C. Chin, J. Hecker Denschlag, and R. Grimm, Sci- 
ence 302, 2101 (2003). 

[10] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. 
Raupach, S. Gupta, Z. Hadzibabic, and W. Ketterle, 
Phys. Rev. Lett. 91, 250401 (2003). 

[11] T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. 
Chevy, M. Teichmann, L. Tarreell, S. J. J. M. F. Kokkel- 
mans, and C. Salomon, Phys. Rev. Lett. 93, 050401 

(2004) . 

[12] C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 

92, 40403 (2004). 
[13] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. 

Raupach, A. J. Kerman, and W. Ketterle, Phys. Rev. 

Lett. 92, 120403 (2004). 
[14] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H. 

Schunck, and W. Ketterle, Nature 435, 1047 (2005). 
[15] D. Jaksch, C. Bruder, J. Cirac, C. Gardiner, and P. 

Zoller, Phys. Rev. Lett. 81, 3108 (1998). 
[16] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Foiling, 

I. Cirac, G. Shlyapnikov, T. Hansen, I. Bloch, and N. 
Contact, Nature 429, 277 (2004). 

[17] T. Kinoshita, T. Wenger, and D. Weiss, Science 305, 
1125 (2004). 

[18] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 



80, 885 (2008). 
[19] D. S. Petrov, Phys. Rev. Lett. 93, 143201 (2004). 
[20] D. S. Petrov, Phys. Rev. A 67, 010703 (2003). 
[21] D. S. Petrov, C. Salomon, and G. V. Shlyapnikov, Phys. 

Rev. Lett. 93, 90404 (2004). 
[22] A. Deltuva and A. Fonseca, Phys. Rev. Lett. 98, 16 

(2007) . 

[23] J. von Stecher and C. H. Greene, Phys. Rev. Lett. 99, 
090402 (2007). 

[24] J. P. D'Incao, S. T. Rittenhouse, N. P. Mehta, and C. H. 

Greene, Phys. Rev. A 79, 030501 (2009). 
[25] Y. Wang and B. Esry, Phys. Rev. Lett. 102, 133201 

(2008) . 

[26] L. D. Faddeev, Zh. Eksp. Teor. Fiz 39, 1459 (1960). 
[27] Y. Suzuki and K. Varga, Stochastic Variational Approach 

to Quantum-Mechanical Few-Body Problems (Springer- 

Verlag, Berlin, 1998). 
[28] R. A. Malfliet and J. A. Tjon, Nucl. Phys. A 127, 161 

(1969). 

[29] O. A. Yakubovsky, Yad. Fiz. 5, 1312 (1967 [Sov. J. Nucl. 

Phys. 5 937 (1967)]). 
[30] J. Macek, J. Phys. B 1, 831 (1968). 

[31] K. Singer, Proc. R. Soc. London, Ser. A 258, 412 (1960). 
[32] S. F. Boys, Proc. R. Soc. London, Ser. A 258, 402 (I960). 
[33] V. I. Kukulin and V. Krasnopol'sky, J. Phys. G 3, 795 
(1977). 

[34] K. Varga and Y. Suzuki, Phys. Rev. A 53, 1907 (1996). 
[35] K. Varga and Y. Suzuki, Phys. Rev. C 52, 2885 (1995). 
[36] K. Varga and Y. Suzuki, Comp. Phys. Comm. 106, 157 
(1997). 

[37] K. Varga, Y. Suzuki, and R. G. Lovas, Nucl. Phys. A 

571, 447 (1994). 
[38] D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev. 

Lett. 99, 233201 (2007). 
[39] J. von Stecher and C. Greene, Phys. Rev. A 75, 22716 

(2007). 

[40] J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. 

A 77, 043619 (2008). 
[41] L. M. Delves, Nucl. Phys. 9, 391 (1959). 
[42] L. M. Delves, Nucl. Phys. 20, 275 (1960). 
[43] C. D. Lin, Phys. Rep. 257, 1 (1995). 
[44] Y. F. Smirnov and K. V. Shitikova, Sov. J. Part. Nucl. 



13 



8, 44 (1977). 

X. Chapuisat, Phys. Rev. A 45, 4277 (1992). [65] 
M. Aymar, C. H. Greene, and E. Luc-Koenig, Rev. Mod. 
Phys. 68, 1015 (1996). 

0. Zatsarinny and K. Bartschat, J. Phys. B 37, 2173 [66 
(2004). [67 

A. Sunderland, C. Noble, V. Burke, and P. Burke, Comp. 
Phys. Commun. 145, 311 (2002). [68 
R. T. Pack and G. A. Parker, J. Chem. Phys. 87, 3888 
(1987). [69 
Y. Zhou, C. D. Lin, and J. Shertzer, J. Phys. B 26, 3937 
(1993). [70 

B. D. Esry, C. D. Lin, and C. H. Greene, Phys. Rev. A [71 
54, 394 (1996). 

H. Suno, B. D. Esry, C. H. Greene, and J. P. Burke Jr, [72 
Phys. Rev. A 65, 42725 (2002). 

F. Werner and Y. Castin, Phys. Rev. A 74, 053604 
(2006). [73 
S. Y. Chang and G. F. Bertsch, Phys. Rev. A 76, 021603 
(R) (2007). 

J. von Stecher, C. H. Greene, and D. Blume, Phys. Rev. [74 
A 76, 053613 (2007). 

Y. Alhassid, G. F. Bertsch, and L. Fang, Phys. Rev. Lett. [75 
100, 230401 (2008). 

1. Stetcu, B. R. Barrett, U. van Kolck, and J. P. Vary, [76 
Phys. Rev. A 76, 063613 (2007). 

H. H. S0rensen, D. V. Fedorov, and A. S. Jensen, Work- 
shop on Nuclei and Mesoscopic Physics: WNMP 2004. [77 
AIP Conference Proceedings, 777, 12 (2005). 
H. H. S0rensen, Arxiv preprint |cond-m at/0502126 (Mas- 
ter Thesis) (2005). [78 
U. Fano, Phys. Rev. A 24, 2402 (1981). 

Even though fi is arbitrary, there is one preferable se- [79 
lection: /i = [Y[ i rm/Q2i mi)] 1/ '' JV_1 ' that conserves the 
volume element of the coordinate transformation. 42] In [80 
this article, we will only present results for equal mass 
systems for which is simpler to simply choose jj, = m. [81 
H. T. Coelho and J. E. Hornos, Phys. Rev. A 43, 6379 
(1991). [82 
N. P. Mehta, S. T. Rittenhouse, J. P. D'Incao, and C. H. 
Greene, Arxiv preprint larXiv:0706.1296l (2007). 
The (3N-7) numerical integration results from the follow- [83 
ing reasoning: initially we have 3N numerical integration 
but 3 dimensions are removed by decoupling the center 
of mass motion, 3 dimension are removed fixing the Euler 



angles and 1 dimension is removed by fixing R. 

J. von Stecher, Ph.D. thesis, University of Colorado, 2008 

available at http://jilawww.colorado.edu/pubs/thesis/ 

vonstecher/von_stecher _thesis.pdf. 

J. Macek, Z. Phys. D 3, 31 (1986). 

D. Blume and C. H. Greene, J. Chem. Phys. 112, 8053 
(2000). 

J. von Stecher, J. P. D'Incao, and C. H. Greene, eprint 

arXiv: 0810.3876 (2008), Nat. Phys. (to appear). 

J. P. D'Incao, J. von Stecher, and C. H. Greene, Arxiv 

preprint larXiv:0903.3348l (2009). 

S. Tan, Arxiv preprint [cond-mat /0412764| (2004) . 

The coefficients p v can be easily related with the s v of 

Ref. [H| by s u =p„_i - 1/2. 

N. P. Mehta, S. T. Rittenhouse, J. P. D'Incao, 
J. von Stecher, and C. H. Greene, Arxiv preprint 
larXiv:0903.4T45l (2009) . 

F. Ferlaino, S. Knoop, M. Berninger, W. Harm, H.-C. 
D'Incao, J. P. Nagerl, and R. Grimm, Phys. Rev. Lett. 
102, 140401 (2009). 

S. B. Papp, J. M. Pino, and C. E. Wieman, Phys. Rev. 
Lett. 101, 040402 (2008). 

C. J. Myatt, E. A. Burt, R. W. Christ, E. A. Cornell, 
and C. E. Wieman, Phys. Rev. Lett. 78, 586 (1997). 

G. Thalhammer, G. Barontini, L. D. Sarlo, J. Catani, F. 
Minardi, and M. Inguscio, Phys. Rev. Lett. 100, 210402 
(2008). 

C. Weber, G. Barontini, J. Catani, G. Thalhammer, M. 
Inguscio, and F. Minardi, Phys. Rev. A 78, 061601(R) 
(2008). 

G. Modugno, G. Roati, F. Riboli, F. Ferlaino, R. Brecha, 

and M. Inguscio, Science 297, 2240 (2002). 

S. Inouye, J. Goldwin, M. L. Olsen, C. Ticknor, J. L. 

Bohn, and D. S. Jin, Phys. Rev. Lett. 93, 183201 (2004). 

M. L. Olsen, J. D. Perre ault, T. D . Cumby, and D. S. 

Jin, Arxiv preprint larXiv:0810. 19651 (2008) . 

T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz, 

and S. Jochim, Phys. Rev. Lett. 101, 203202 (2008). 

J. H. Huckans, J. R. Williams, E. L. Hazlett, R. W. Stites, 

and K. M. O'Hara, Arxiv preprint larXiv:0810.3288l 

(2008). 

D. Blume, S. T. Rittenhouse, J. von Stecher, and C. H. 
Greene, Phys. Rev. A 77, 33627 (2008). 



